import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate

# 中文和负号的正常显示
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 壳层法计算旋转体体积
def shell_method(f, a, b, axis='y'):
    """
    使用壳层法计算旋转体体积
    f: 函数
    a, b: 积分区间
    axis: 旋转轴，'y'表示绕y轴旋转，'x'表示绕x轴旋转
    """
    if axis == 'y':
        # 绕y轴旋转，对x积分
        volume, error = integrate.quad(lambda x: 2 * np.pi * x * f(x), a, b)
    else:
        # 绕x轴旋转，对y积分
        # 这里假设f是x关于y的函数
        volume, error = integrate.quad(lambda y: 2 * np.pi * y * f(y), a, b)
    
    return volume, error

# 示例：y = x^2 在 [0,2] 上绕y轴旋转
def f2(x):
    return x**2

volume_shell, error_shell = shell_method(f2, 0, 2, 'y')

# 使用圆盘法验证
def inverse_f2(y):
    return np.sqrt(y)  # x = sqrt(y)

volume_disk, error_disk = integrate.quad(lambda y: np.pi * (inverse_f2(y))**2, 0, 4)

print(f"壳层法计算的体积: {volume_shell:.6f}")
print(f"圆盘法计算的体积: {volume_disk:.6f}")
print(f"理论值: {8*np.pi:.6f}")

# 可视化壳层法原理
x = np.linspace(0.1, 2, 20)  # 避免x=0
theta = np.linspace(0, 2*np.pi, 50)

fig = plt.figure(figsize=(10, 5))
ax = fig.add_subplot(111, projection='3d')

# 绘制几个代表性的圆柱壳
for i, xi in enumerate(x[::4]):  # 每隔4个取一个点
    theta_i = np.linspace(0, 2*np.pi, 30)
    y_i = f2(xi) * np.ones_like(theta_i)
    x_i = xi * np.cos(theta_i)
    z_i = xi * np.sin(theta_i)
    
    ax.plot(x_i, y_i, z_i, 'b-', alpha=0.7)

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('壳层法原理示意图')

plt.show()